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FOREWORD 


A knowledge of the energy influx along a shoreline is of consider- 
able importance to an understanding of many beach processes and long- 
term trends in beach responses. This study represents a step toward the 
development of a rapid method for routine determinations of wave power 
along a shoreline, using observed or hindcast deep-water wave character- 
istics and high-speed computer programs for the calculation of wave 
refraction. Specifically treated here is a computer program and pro- 
cedure for the construction of wave rays. An example of the method is 
presented in which wave rays are brought from deep water into the Atlantic 
shoreline of the City of Virginia Beach, Virginia. Detailed explanations 
of the computer programs used, instructions for their operation, and 
sample inputs and outputs are given in appendices. A series of sugges- 
tions is also given for improvement of the method presented. 


This report was prepared by Dr. Wyman Harrison, Associate Marine 
Scientist, Virginia Institute of Marine Science, in pursuance of Contract 
DA-49-055-—CIV-ENG-64-5 with the Coastal Engineering Research Center, and 
in collaboration with Mr. W.Stanley Wilson, formerly a graduate student 
at the Institute. 


The study was supported by the Coastal Engineering Research Center 
(formerly the Beach Erosion Board of the Corps of Engineers). Computing 
Centers at Northwestern University, the College of William and Mary, and 
NASA (Langley Field, Virginia) extended every cooperation. Lieutenant 
G. Griswold, Oregon State University, first suggested to the authors the 
use of numerical methods for the calculation of wave refraction. The 
computer program used in this report for calculating wave refraction was 
extensively modified from a program under development in 1963 by Lieu- 
tenant G. Griswold and Mr. F. Nagle of the U. S. Navy Weather Research 
Facility, Norfolk, Virginia, and Mr. E. Mehr of New York University. 

Mr. J. Gaskin of IBM and LCDR C. Barteau of the U. S. Navy Weather 
Research Facility aided in adapting the Griswold-Mehr program for use on 
the 7094 computer, while Mr. J. Curran and Mr. R. Libutti of IBM, aided 
in attempts to adapt the program to the 1620 computer. Both programs 
were considerably revised by Wilson prior to their presentation in the 
appendices of this study. 

The authors are indebted to Professor W. C. Krumbein of Northwestern 
University, Consultant to the Coastal Engineering Research Center, and to 
Mrs. Betty Benson of the Northwestern University Computing Center, for 
making available the surface-fitting program that is used in the PRMAT 
subprogram and SURFCE subroutine of the wave-refraction program. Pro- 
fessor B. R. Cato of the Mathematical Department of the College of 
William and Mary, provided helpful advice at various stages of the study. 


This report is published under authority of Public Law 166, 79th 
Congress, approved July 31, 1945, as supplemented by Public Law 172, 
88th Congress, approved 7 November 1963. 
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DEVELOPMENT OF A METHOD FOR NUMERICAL 
CALCULATION OF WAVE REFRACTION 


by 


W. Harrison 
U. S. Coast and Geodetic Survey, Washington, D. Ca and 
Virginia Institute of Marine Science 
and 
W. S. Wilson 
The Johns Hopkins University, Baltimore, Maryland? 


ABSTRACT 


Steps in wave-ray construction are as follows: 


1. Select wave periods and approach angles for each 
series of rays to be constructed. 


2. Prepare a grid of depth values for the area of 
investigation. 


3. Use a computer program to obtain a table of 
water depths and related wave velocities for each wave. period 
Sellefetedmin Ge) 


4, Make a grid of wave-velocity values for each 
wave period selected in (1), for the area of investigation, 
using a second computer program which takes as input the depth 
grid of (2) and the appropriate depth-velocity table of (3). 


5. Derive, using a third computer program, matrices 
for use in the linear interpolation scheme of the computer 
program of (6). 

6. Calculate, for each wave period specified in (1), 
the points along a wave ray using a computer program which 
takes for input: 

a. The appropriate velocity grid of (4). 
b. The matrices of (5). 


c. The origin points and approach angles of 
(1) for given wave rays. 


lpresent address. Study completed while at the Virginia 
Institute of Marine Science. 


A linear-interpolation scheme (using the least-squares 
method) is used in determination of wave velocity at a given 
point along a ray. Ray curvature is then calculated at this 
point and an iteration procedure is solved to obtain the 
position of the next point. The ray terminates at the shore 
or grid border. 


An example of usage of the programs is presented in which 
52 rays for waves of 4- and 6-second period, from six different 
deep-water origin points, are brought into the Atlantic shore- 
line of the City of Virginia Beach, Virginia. 


The procedure outlined is in the developmental stage, 
and suggestions for improvements are given that should offer 
a quick, accurate, and objective method of constructing wave 
rays. 


INTRODUCTION 
Background of Project 


An ultimate understanding of the changes in topography of a given 
shoreline will be derived from a considerably fuller knowledge of the 
input and transformation of energy along the strand than is currently 
available. Because by far the greatest amount of energy expended in the 
beach-ocean-atmosphere system is associated with breaking waves in the 
surf zone, it becomes of the greatest importance to evaluate the long- 
term distribution of wave energy along the shoreline. This energy dis- 
tribution depends mainly upon the energy of the original deep-water waves, 
as modified by refraction when the waves move shoreward through shallow 
water. 


It was with these general considerations in mind that the authors 
embarked upon a study of wave refraction at Virginia Beach, Virginia, 
where previous studies (cf. U. S. Congress, 1953; Harrison and Wagner, 
1964) suggested the desirability of determining wave power distributions. 
The authors were armed with an “operational"' computer program said to be 
capable of rapid and accurate calculation of wave-ray paths. As in the 
case of many "operational" computer programs (cf. the cogent comments of 
Adams, 1964), the authors soon discovered that they were either in diffi- 
culty simply trying to make the program operate, or that they were not in 
agreement with certain of the logical steps involved. It soon became 
apparent that the entire method needed to be reviewed and revised. 


The emphasis here on development of a numerical method for calcu- 
lation of wave refraction does not mean that the long-range goal of 


assessing wave-power distributions along shorelines has beén set aside. 
The authors feel, moreover, that the high-speed method ultimately adopted 
for calculating wave refraction is of such basic importance to the larger 
problem that a detailed presentation is warranted at this time. 


Earlier Studies of Wave Refraction 


The investigations of Krumbein (1944) and Munk and Traylor (1947) 
were among the first in which steps were taken to assess the link between 
wave refraction, the energy distribution along the shore, and beach erosion 
or deposition. These workers constructed wave-refraction diagrams for 
selected deep-water wave periods by hand-drawn methods, a time-consuming 
process that at best permits only partial assessment of the effect of a 
spectrum of waves on shore erosion. 


About ten years ago, Pierson, Neumann, and James (1953) considered 
wave refraction effects in some detail and concluded (1953, p. 186) that 
to use only significant height and the average "period" of the waves in 
deep water and to refract the waves with these two numbers will lead to 
totally unrealistic results. These authors went on to outline a method 
(1953, p. 197) for forecasting wave heights and characteristics at a point 
in shallow water near a coast. The method involves construction of a 
‘large number of ray diagrams by graphical methods. 


Except for the great amount of work involved, the method holds 
promise for a realistic assessment of energy distributions along a shore- 
line, assuming the deep-water spectra can be adequately approximated. 
Although more "practical" methods for determining wave refraction (cf. 
Silvester, 1963) and for assessing wave energy along coasts (cf. Walsh, 
Reid, and Bader, 1962) have been proposed, the present authors feel that 
the adequate treatment of the problem of energy inflow must of necessity 
include large-scale computations. The treatment of wave refraction, for 
example, should be as thorough as that outlined by Pierson, Neumann, and 
James (1952, p. 197), and this requires many man hours. The advent of 
high-speed computers offers the possibility of significant reductions in 
the time required for construction of refraction diagrams, wave spectral 
forecasts, and, resultant distributions of energy at the shore. The 
machine methods given below for constructing refraction diagrams represent 
a step in the attempt to develop a comprehensive computer program for the 
routine determination of wave energy along shorelines. 


Scope of Present Study 
It will be assumed that the deep-water wave spectral periods, heights, 


and directions are known, and that they may serve as input to a computer 
program for determining wave refraction as part of an overall method such 


as proposed by Pierson, Neumann, and James (1953). Examples of the 
machine refraction of six separate combinations of wave period and 
direction off Virginia Beach, Virginia, are given to illustrate the 
method of refracting wave rays. The assumptions and methodology used 

in constructing the wave rays are fully treated, and suggestions for im- 
proving and amplifying the program are then presented. 


METHOD 
Previous Work 


The presentiy acceptable method for constructing wave-ray diagrams 
involves the hand method of graphical construction first proposed by 
Arthur, Munk, and Isaacs (1952) and later repeated in Pierson, Neumann, 
and James (1953) and U. S. Army, Corps of Engineers (1961) Technical 
Report No. 4, of the Beach Erosion Board. This procedure for graphical 
construction of wave rays involves the use of Snell's Law, 


sin Oy 


sin oP) Cy 


where for two successive points ona ray, 1 and 2, C = wave velocity, and 
a = angle between the wave crest and the bottom contours. In the hand- 
refracting methods, Snell's Law requires the assumption that between points 
1 and 2 there exist straight and parallel contours. The computer program 
for-construction of wave rays (that is to be described later) involves the 
fitting of interpolative surfaces to points of a wave-velocity grid. This 
program, although requiring certain assumptions of its own, avoids the 
assumption of straight and parallel contours between successive points on 

a ray. With this program an expression for ray curvature is solved quite 
independently of the use of Snell's Law. 


The authors' introduction to the possibilities of computer programs 
for the numerical construction of wave rays was gained originally through 
discussions with Lieutenant Gale M. Griswold, then of the U. S. Navy 
Weather Research Facility, Norfolk, Virginia. Ideas on the subject had 
appeared in mimeographed reports (Mehr, 1961, 1962a, 1962b; Griswold and 
Nagle, 1962) and a published paper (Griswold, LOGS). A computer program, 
although not actually operational, was provided by Griswold for develop- 
ment in this study. (This program will be referred to henceforth as the 
Griswold-Mehr program. ) 


Construction of Wave Rays 


Selection of Input. The first step in wave-ray determination begins 
with the construction of a grid of depth values for the area of interest. 


Geographical limits for a grid of depth values take into account the 
probable origin points for waves of all deep-water (d/Lo>>0.5) approach 
angles of interest. By convention, the grid origin is located in the 
southwestern corner of the area to be covered, with the X axis extending 
eastward along the southern border and the Y axis extending northward 
along the western border. The grid interval is selected such that, in a 
given cell, bottom contours can be represented by straight and parallel 
lines. (This representation will be discussed in detail later.) 


Actual or interpolated depths at grid intersections are recorded to 
the nearest foot, and all real depths are made positive. Contours are 
then drawn of depth values extending several grid units out from shore. 

At this point, contours that are symmetrical reflections of nearshore 
bottom contours are drawn on the land, over a perpendicular distance from 
the shore of two grid units. "Depths" at the grid points on land are then 
derived in the region of the symmetry contours on land; these "depth" 
values are made negative. All other land "depths" (that is, those farther 
than two grid units from the shore) are made zero. (The procedure of 
assigning negative values for nearshore "depths" on land is required for 
the fitting of wave-velocity surfaces.) The depth values so obtained are 
prepared as input to the DISTV program (Appendix B). 


Example of Input. An area of the mid-Atlantic Bight (Figure 1) is 
chosen to illustrate the method of data input and the computation of wave 
refraction. A specific wave-ray target area if the region covered by the 
depth grid is the shoreline of the Borough of Virginia Beach in the City 
of Virginia Beach, Virginia (Figure 1). 


The wave input data used here are not representative of any given 
wave spectrum observed or hindcast. They approximate the 15 most commonly 
observed combinations of wave height, period, and direction observed at 
the Chesapeake Lightship (Figure 1). These combinations were determined 
and condensed for input to the computer at the request of the Coastal 
Engineering Research Center. The method used in determination of the 
combinations is described in Appendix G. The authors have used the six 
wave period -- direction combinations (Tables 1 and 2), extracted from 
the 15 height-period-direction combinations of Appendix G, merely to 
illustrate the method of ray construction. 


The depth grid (Figure 1) of the example was chosen so that the 
origin points of six second waves approaching Virginia Beach from 60° and 
90° T would be positioned in water depths of greater than 92 feet 
(d/Lo>0.5 for T = 6 sec.). 


Outlines of the 99 X 81 unit depth grid used in the example are 
shown on Figure 1; the grid origin is located at 7691.9'W, 36°39.5'N, and 
the grid interval is 3,040 feet. U. S. Coast and Geodetic Survey (boat- 
sheet) charts 5988, 5990, 5991, 5992, 5993, and 6595 were used in obtaining 


depth values. Where these charts did not offer coverage, depths were 
picked off charts 1222 and 1227. Smoothed contours of the depth values 
at grid points appear on Figure 2. 


Deep-water starting points and directions for the 52 wave rays of 
the example are given in Tables 1 and 2, and are shown plotted on Figure 
3. An example of the coding of thewave-ray and wave-velocity input data 
is given under the heading "INPUT" in Appendix D. 


Computer Operations 


Programs Used. For each wave period selected for the investigation, 
a table of depths and their associated wave velocities is prepared. This 
is accomplished by solving the following equation (U. S. Army, Corps of 
Engineers, 1942; U. S. Navy, Hydrographic Office, 1944; Mason, 1950) by 
an iteration process: 


_ 2d 
C= TT tanh Tic) 


where C = wave velocity, g = acceleration due to gravity, T = wave period, 
and d = water depth. This equation has been programmed, and a sample 
depth-velocity table has been prepared for T = 4 seconds. (See COMPV in 
Appendix A.) 


As in other wave refraction studies (Pocinki, 1950; Pierson, 1951; 
Pierson, Neumann, and James, 1953), it is assumed that wave velocity is a 
function only of water depth and wave period, as expressed by the above 
equation. Various factors such as bottom friction, percolation, reflection, 
currents, and winds are considered as having no effect on the refracting 
waves. 


Given the absolute value of a water depth, it is possible to check 
the appropriate table, which has been previously prepared, for the 
associated wave velocity. Preparation of an entire velocity grid for 
each wave period to be studied is then carried out. This procedure has 
also been programmed, and a portion (from X =.0 to X = 19, and from 
Y = 0 to Y = 2) of the depth grid described above has been derived for 
T = 4 seconds. (See DISTV in Appendix B.) 


The procedure for determining the wave-velocity value at an 
arbitrary point in a grid cell involves fitting a plane to the velocity 
values at the four grid intersections that bound the grid cell in question. 
This linear surface is fit by the least-squares method, using an equation 
of the form: 


Cr Et BOX eB, 


where C = wave velocity, E's = coefficients, and X,Y are the grid coor- 
dinates for the arbitrary point. 


With the least-squares method of surface fitting, it is possible to 
obtain certain matrices which are used each time the equation of a plane 
is derived for a grid cell. These matrices have been derived in PRMAT. 
(See Appendix C.) That portion of the surface-fitting procedure which must 
be carried out each time a linear equation is to be derived is given in 
SURFCE subroutine of MAIN 1620 (Appendix D) and MAIN 7094 (Appendix E). 


Determination of each desired ray for a given period is accomplished 
by first specifying origin coordinates and an angle of approach. The actual 
ray is constructed by plotting and connecting the series of successive 
calculation points (computed by MAIN 1620, Appendix D, or MAIN 7094, Appendix 
E) which range across the velocity grid until striking the beach or grid 
margin. (As discussed by Griswold (1963, p. 1722) the rays may be run from 
the shore seaward, if the construction of a refraction diagram at a point 
is desired.) Velocity at each point is calculated as mentioned above; ray 
curvature (FK) is calculated by using the following expression (Munk and 


Arthur, 1952): 
FK = = | sin A (&) - cos A () | 


where A = approach angle, as defined in Appendix D. 


In order to determine Xpn+1, Yn+i4, An+1, and FKyn+ 1 for calculation point 
n+1 (with those similar values known for point n), the following equations 
(Griswold and Nagle, 1962; Griswold, 1963) are solved by an iteration 
procedure: 


AA = (FE + FK 40/2 
A =A + AA 

n+1 n 
A = (AL + Aled /2 
Xie =X + DcosA 
vee = Y +DsinA 


where D equals the incremented distance between points n and n+l. (See 
MOVE subroutine of MAIN 1620, Appendix D, and of MAIN 7094, Appendix E.) 


The computers to be used for the programs cited in the appendices 
are the IBM 1620 (with floating-point hardware and 60K memory) and the 
IBM 7094. FORTRAN II is the language used for the 1620 programs; FORTRAN 
Iv is used for the 7094 programs. The computer and language for a given 
program are specified on the first comment card of each program print-out. 


The number of digits carried by a computer during internal computations 
for floating-point numbers (designated by f) and the number of digits carried 
for fixed-point numbers (designated by k) vary with the different programs. 
For the 1620 programs, f is given in Columns 2 and 3 and k is given in 
Columns 4 and 5 of the first card of each program print-out. For the 7094 
program, f = 8 and k = 5. 


Example of Sample Output. Sample input and output listings are given 
in Appendix D for three wave rays. These three rays, numbered 1, 2, and 3 
in the OUTPUT listings, correspond to rays numbered 33, 32, and 31, respec-— 
tively, of Table 1 and Figures 2 and 6. 


Representation of Output 


Possible Methods. Output from the computer programs could be presented 
in a variety of ways, depending upon the requirements of the investigator. 
Precision work on caustics, for example, might entail skillful plotting of 
the values at MAX points with precision drafting techniques. If a rapid 
analysis of ray paths from deep water to the shore is all that is required, 
however, the investigator might consider an X-Y plotter attachment for 
rapid printing of the computer output. For many shore engineering studies, 
only the final few hundred yards of the wave rays will be of practical 
value, while for studies of energy or wave-power distribution along a 
fixed segment of the shore, it is possible that only the terminal points 
of the rays at some specified distance from shore would be of value for 
presentation. 


Example. The results of the computer refraction of the 52 wave rays 
of Figure 2 are presented in Figures 3-8, which show only the last portions 
of the wave rays relatively close to shore. Successive calculation points 
were plotted on a grid and then connected by a smooth curve. The rays 
themselves vary from the almost unmodified ones’for 4-second waves that 
approach the beach relatively head on (Figure 5), to the 6-second ones 
that actually cross (Figure 8) within the limits of the figure. 


Attention is drawn to the terminal points of the wave rays; some of 
these points (such as those of rays Nos. 14, 37, and 52 as shown in Figures 
4, 7, and 8, respectively) are closer to the shoreline than others. These 
variations are due to the fact that a ray terminates because either the 
next calculation point is in an area of zero or negative velocity, or 
curvature approximations are not converging. Thus, with a constant in- 
cremented distance (D), all rays do not reach a similar distance from the 
shoreline. Ray No. 24 (Figure 5), on the other hand, makes a pronounced 


curve near its terminus that appears out of context when compared to the 

adjacent wave rays. This is due to a combination of (1) a poor "fit" by 

the linear-interpolation surface at this point, and (2) poor grid control 
(i.e., poor representation of depth values recorded on the initial depth 

grid). Suggestions for elimination of these factors are discussed in the 
section "Grid Considerations." 


PROGRAM DEVELOPMENT 


The computer programs presented in the appendices have by no means 
been refined to the fullest extent. At the moment the following factors 
may be considered of most importance to their continued development: (1) 
improvement of the surface-fitting scheme for wave velocity interpolation, 
(2) improvement of the method of ray-curvature approximation, (3) addition 
of a provision for changes in grid scale and incremented distance as a ray 
moves shoreward, (4) refinement of ray-plotting techniques, and (5) testing 
of the revised program at a suitable place in nature. 


Interpolation Surfaces 


"Forced-cubic" Interpolation Surface. First mentioned by way of 
review is the surface-fitting scheme for interpolation of wave velocity 
used in the Griswold-Mehr program. Their scheme involves fitting a cubic 
surface that (1) passes exactly through the velocity values at the grid 
intersections of the given grid cell, and (2) is the cubic surface of 
best fit (by the least-squares method) to the velocity values at the 
eight grid intersections closest to and surrounding the four grid inter- 
sections of the given cell. This cubic surface is called a "forced-—cubic" 
surface because it is "forced" to pass through the four innermost velocity 
values. Because it permits the taking of first and second derivatives for 
use in wave intensity calculations, as explained by Griswold (1963, p. 1720), 
it was the first surface-fitting program used in this study. 


An example of the results of its use appears in Figure 9C. It is 
obvious from the figure that this method of interpolation is invalid. The 
forced feature of the surface creates undesired ridges and/or troughs in 
the velocity surface. Thus, depending upon the location of a given calcu- 
lation point in a grid cell, the ray can be deflected erratically from a 
"normal" path. Results such as those exemplified in Figure 9C necessitated 
a search for an interpolation surface that could adequately portray the 
general trend of wave velocity change in a given cell. 


Quadratic-interpolation Surface. The Griswold-Mehr program was 
altered by insertion of a quadratic surface of best fit (by the least- 
Squares method) subroutine for preparation of the interpolation surface. 
In order to derive a quadratic surface, at least six data points are 
necessary. It was decided to use the velocity values at the closest 12 
grid intersections. The results (Figure 9B), however, did not yield a 


satisfactory interpolation surface; this was evident when calculation 
points for a ray had moved within two grid units of the shore. This is 
because the negative "land" velocity values, used in derivation of the 
surface, made the general tilt of the surface steeper than indicated by 

the velocity values at the central four grid intersections. This is 

shown in Figure 9B where the rays bend, after moving within two grid 

units of the shore, more than do the rays produced by the linear-interpola- 
tion program (Figure 9A) or the rays produced by the graphical-construction 
method (Figure 9D). 


Linear-interpolation Surface. The linear-interpolation programs 
(MAIN 1620, Appendix D; and MAIN 7094, Appendix E) were developed in order 
to remedy the excessive tilt created by the quadratic-interpolation program. 
The advantage offered by the linear-interpolation program is that only the 
velocity values at the central four grid intersections are needed in order 
to derive the surface. The assumption, presented under the heading "Selec- 
tion of Input," stating that the grid interval be selected such that bottom 
contours in any grid cell be represented by straight and parallel lines, is 
not actually valid. Although bottom topography may be represented by a 
plane (i.e., it changes linearly) in a given cell, the associated velocity 
values may be changing exponentially (i.e., velocity is an exponential 
function of depth) in the same cell. However, for purposes of this program, 
it is assumed that the velocity values in a given grid cell can be adequately 
represented by a plane. 


The use of the variable PCTDIF in MAIN 1620 (card number RAYN 19, 
Appendix D) and MAIN 7094 (card number RAYN 20, Appendix E) serves to give 
an indication of the relative "fit" of a surface to a given set of velocity 
values. The output values for PCTDIF (see OUTPUT, Appendix D) give the 
percent difference between only one of the four velocity values at the 
nearest four grid intersections and its related value on the plane fit to 
the same four values. PCTDIF by no means represents the degree of "fit" 
of the surface to the four values because these percent differences may 
vary for each of the four values to which a plane is fit. This is especially 
true in cells where there are great differences between the four velocity 
values (such as near the shore). PCIDIF does, however, give an estimate of 
the relative error encountered in the interpolation in a given grid cell. 


Just as the quadratic-interpolation program yielded an excessive tilt 
when the given grid cell in use fell within two grid units of the shore, 
the linear-interpolation program yields an excessive tilt when the given 
cell falls within one grid unit of the shore. Therefore, it is expected 
that PCTDIF will be large in value when the grid cell being interpolated 
is located near the shore. In view of this fact, the rays run with the 
linear surface-fitting program should be considered as rough approximations 
only, when closer than one grid unit to the shore. 


Interpolation Surface Versus Graphical Method of Ray Construction. 
Consideration of the above-mentioned surface-fitting programs led to the 
choice of the linear surface of best fit as the most suitable one for use 
in construction of wave rays. In general, it was found that rays run with 
the linear surface of best fit compared most favorably with rays constructed 
by graphical methods, and this is the case in the plotting example (Figure 9 
A and D). An advantage of the least-squares linear-interpolation method 
over the graphical method lies in the fact that ray curvature can be com- 
puted at a number of points between contours. Aside from the absolute 
validity of the two methods of wave-ray construction, the computer method 
is estimated to be many times faster than the hand method when only the 
terminal points are desired of a large number of rays. Where only a few 
rays are desired, the hand method is clearly the most rapid and economical. 
Because the practice of refracting only a few rays for a few wave periods 
yields totally inadequate information upon which to assess wave heights and 
energies at the shore, it seems clear that the real considerations are not 
so much the man-hours involved in depth-grid and program-deck preparation 
for the computer as they are in the necessity for the more realistic results 
that can be afforded by the computer construction of wave rays for entire 
wave spectra. 
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Future Interpolation Schemes. A more valid interpolation scheme 

would be obtained if a grid of depth values were input to MAIN 1620 or 
MAIN 7094 instead of a grid of velocity values. The assumption, that the 
grid interval be selected such that depth contours in any grid cell be 
represented by straight and parallel contours, would then be perfectly 
valid. In this case, a depth value would be obtained from the interpola- 
tion surface at a given calculation point. Then, using the procedure used 
by COMPV (Appendix A), the depth value could be converted into a wave- 
velocity value. It is noted that, in order to obtain curvature (FK) at a 
calculation point, && and 35 are needed. If a linear-interpolation program 


Ge used, with an input depth grid, expressions would be available for 


x and a: on oor ou tne. ne Lact onsntp (derived in Appendix F) could then 
e used to obtain 2+ andl: 
x oy 
oc = od ZA le = od Z 
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In this expression k' = T/4T and k" = 21/gT. If a quadratic surface were 
to be used for interpolation similar relations could be derived to relate 
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A reason for using a quadratic (or a higher-order polynomial) surface 
is that the second partial derivatives can be obtained. This fact was 
used in the Griswold-Mehr program to obtain values for Beta (Munk and 
Arthur, 1952), the coefficient of refraction, and H/Ho. (As mentioned, 
however, their interpolation scheme rendered their results invalid.) It 
is apparent that certain problems with the quadratic (ar higher-order 
polynomial) interpolation schemes must be resolved before such sophisticated 
parameters as Beta or H/Ho can be estimated. It should be possible to 
derive a quadratic surface by heavily weighting the velocity values at the 
central four grid intersections while still using the surrounding eight 
velocity values. This procedure may produce a good interpolation surface; 
if so, it would be quite worth while to calculate estimates of the above 
parameters. 


Ray Curvature Approximation 


An entire grid of velocity values is not a smooth and continuous 
surface when observed as a set of planes in which each grid cell is 
represented by a specific plane of best fit to its four bounding velocity 
values. With this in mind, it is not surprising that all curvature approxi- 
mations (see MOVE subroutine of MAIN 1620, Appendix D; or MOVE subroutine of 
MAIN 7094, Appendix E) fail to converge to a single value when determining 
a new point along a ray path. This is especially apparent when adjacent 
planes present a large discontinuity in wave-velocity values at their 
common edges. This fact is the reason for the variable MIT included in 
MOVE subroutine. In case the curvature approximations oscillate among 
three or more values after 20 iterations, the ray is terminated, as no 
valid curvature approximation can be made. Although infrequent, sometimes 
(not shown in OUTPUT, Appendix D) the curvature approximations oscillate 
between two values. In this case, the message "CURVATURE APPROXIMATED" is 
included with the output in order that the operator note that the curvature 
used for the given point was an average of two values. The Griswold-Mehr 
program did not have such a check on the curvature approximations; the 
average of the last two approximations after 10 iterations was used as the 
new curvature value, if the values did not converge. This is another reason 
for the erratic ray behavior near the shore in Figure 9C. 


Grid Considerations 


For reasons outlined in a previous section ("Future Interpolation 
Schemes") assume that a grid of depth values will be input and used for 
derivation of the interpolation surfaces. It will be necessary in the 
future program to transfer a ray to a larger scale grid when approaching 
the shore in order to allow better grid control (i.e., better representation 
of depth values) in an area where the velocity values are rapidly changing. 
There is, however, a limit to the maximal scale of a grid as beach and 
nearshore topography are constantly changing, especially during each storm. 
Thus the selection of grid interval is a question to be pursued in additional 
studies. 


It will also be necessary to vary the distance (D) incremented 
between successive calculation points along a ray in order that the ray 
may proach as close to the shoreline as possible. As an example, the 
value of D could be assigned such that D = 0.5 when d/Lg>0.5 and D = d/Ly 
when d/L .<0.5. 


Wave Ray Representation 


At present, output must be plotted by hand; and, if there are a 
large number of rays to be constructed, this can be a tedious and time- 
consuming process. The use of an X-Y plotter, if adaptable to changing 
grid scales, would be an ideal scheme for rapid plotting of wave rays. 
Another scheme might involve the hand or machine plotting of the (numbered) 
terminal points of the wave rays, along the shore or grid margin. 


Testing of Ray Constructions 


The absolute validity of the ray constructions (Figures 3-8) is 
impossible to evaluate without a test in nature. It would be desirable 
to test the ray constructions in an area where significant refraction of 
relatively constant-period wave trains occur and where the bathymetry is 
well known. The use of aerial photography would be an invaluable aid in 
conducting such a test. 


It is possible to test ray constructions by comparison with analytic 
solutions for wave rays passing over algebraically-described surfaces 
(Pocinki, 1950). Griswold (1953), in testing the Griswold-Mehr program, 
used both a straight uniformly-sloping beach and a parabolic bay and found 
little variation between computed and analytic rays. However, such a 
theoretical test is not sufficient reason for acceptance of a computer 
program utilizing an interpolation scheme. In these theoretical tests, 
an algebraically-described surface of velocity values is input to the 
computer. The close agreement then noted between computed and analytic 
rays is due to the fact that the given interpolation scheme can easily and 
accurately represent the velocity surface when interpolating within any 
grid cell. In order to adequately test such a computer program, an 
irregular surface of velocity values must be provided. This necessity 
supports the need for a test in nature, as mentioned in the previous 
paragraph. 


CONCLUSION 


The procedure described in this report, when fully improved using 
the accompanying suggestions, will constitute a rapid and accurate method 
of wave ray construction. At the present stage of development, however, 
the procedure must be accepted with reservation. 
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TABLE 1.-Computer input specifications for waves of 4~sec period 


Ray Ray origin Grid origin Direction from which 
no. angle for ray travels 
computer 
(A) Cos) (Cr) 
i 20.0 Bo 50) Go 50 30 
2 20.0 RsAs (O00 30 
3 20.0 39.25 75.50 30 
4 240.0 52 7/5100 30 
5 2h0 .0 LOO Who 5O 30 
6 20.0 MLa i FAC 30 
i 20.0 ne) ih Ges} 0) 30 
8 210.0 2659239. 23 60 
9 210.0 PMS 3S38 60 
10 210.0 BO Si 52 60 
afin 210.0 28.45 36.66 60 
12 210.0 28.96 35.80 60 
ali 210.0 29.47 34.94 60 
14 210.0 29.98 34.08 60 
15 210.0 BOgMe) 3-22 60 
16 210.0 31.00) 32.36 60 
iY 210.0 SL G0) Subs SO 60 
18 180.0 250 26.150 90 
19 180.0 PAO 25550 90 
20 180.0 250) 2.50 90 
Pill 180.0 21.50 23.50 90 
22 180.0 Pilie50 22,50 90 
23 180.0 PAL 50) 2, SO) 90 
pn 180.0 PS Om ee2 On © 90 
25 180.0 P50). 1G. 50 90 
26 180.0 250 UIo50 90 
oT 180.0 Pio 50) Io 50 90 
28 120.0 siscil OHC6 150 
29 120.0 oO Oi 150 
30 120.0 Oks) O84 Si/ 150 
31 120.0 1622 OSES 150 
32 120.0 15.36 02.99 150 
33 120.0 14.50 02.50 150 


Table 2.-Computer input specifications for waves of 6-sec period. 


Ray Ray origin Grid origin Direction from which 
no. angle for ray travels 
computer 
(A) (x) (x) (ME) 

34 210.0 85.50 73.50 60 
35 210.0 SSO, 7Aseh 60 
36 210.0 G52 Wil 60 
37 210.0 87.03 70.92 60 
38 210.0 87.54 70.06 60 
39 210.0 88.05 69.20 60 
ho 210.0 88.56 68.34 60 
Ta 210.0 89.07 67.48 60 
he 210.0 89.58 66.62 60 
43 180.0 Ci, 50  2S.50 90 
yh SOO Ci, 50 25050 90 
45 180.0 CisSO) Des 50 90 
h6 180.0 91.50 23.50 90 
LT 180.0 91.50 22.50 90 
48 180.0 GiLe5O ilo 50 90 
hg 180.0 91.50 20.50 90 
50 180.0 91.50 19.50 90 
bil 180.0 Co 50) Msig 59 90 
52 180.0 91.50 17.50 90 
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FIGURES 
Text 
Map of portion of mid-Atlantic bight showing generalized bathymetry and 
areas covered by depth grid and detailed ray diagrams. 
Map showing contours of depth values associated with intersection points 
of the primary grid, and the starting points and directions of the 52 
wave rays run with the computer programs. 
Wave-ray diagram for rays numbered 1-7 (fig. 2). 
Wave-ray diagram for rays numbered 8-17 (fig. 2). 
Wave-ray diagram for rays numbered 18-27 (fig. 2). 
Wave-ray diagram for rays numbered 28-33 (fig. 2). 
Wave-ray diagram for rays numbered 34-2 (fig. 2). 
Wave-ray diagram for rays numbered 43-52 (fig. 2). 
Comparisons of wave rays 31-33 (fig. 2) constructed by three different 
programs for fitting velocity surfaces (the linear-surface program, A; 
the quadratic-surface program, B; and the forced-cubic-surface program, 
C) and by conventional graphical methods, D. 
Appendix G 
Map showing former location of U. S. Navy wave gage off Cape Henry, 


Virginia, in 20 feet of water. 
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APPENDIX A 
Computer Program for Computation of Wave Velocities as a 
Function of Wave Period and Water Depth 
PROGRAM TITLE: COMPV. 
Variables Used in Program: 
INCHES So d5000 --- Number of different wave periods for which velocity computa- 
tions are to be made. 
Mente eae npeciodm (seconds) 
KADEP eerie Wateridep thm cnecty re 
CXX(K)........ Wave velocity (feet/second) for water of K depth. 
XL.......2.2-- 1/2 deep-water wave length (feet) for a wave of TT period. 
L..eeeee0e-ee-- XL YOUNnded down to the nearest integer. 
Summary of Program: 
NOTT and the first TT are the first variables input to the computer. For 
water of K depth, where K ranges from 1 to L in one-foot increments, CXX's 
are then computed by an iteration pro¢edure so that the third digit to the 
right of the decimal is significant. After all CXX's have been computed 
for the first TT, they are rounded to the nearest hundredth and output. 
Computations then proceed using the next TT. After computations have been 
made using NOTT different TT's, the program stops. 
Remarks: 
If wave periods greater than 20 seconds are to be used with this program, 
CXX will require larger dimensions. 
The output from this program consists of a table of water depths and asso 
ciated wave velocities for each specified wave period. The CXX values of 


these serve as input for the DISTV program in Appendix B. 


30 


The COMPV program, as well as several of the following programs, utilizes a 
subroutine for the internal rounding of values before they are output. This 
subroutine (ROUND) may profitably be described at this point. 
SUBROUTINE TITLE: ROUND .« 
Variables Used in Subroutine: 
VALUE.......-. Corresponds to the value in the calling program which is to 
be output. 
WIRGooooGCOK0R0 10% where N is the number of digits which are to be output 
to the right of the decimal. 
Summary of Subroutine: 
If there are N digits to the right of the decimal in the output specification 
for VALUE, ROUND tests the N+lst digit to see if it is equal to or greater 
than 5. If it is, the Nth digit is rounded up. Otherwise, the Nth digit 
remains unchanged. | 
Remarks: 
Mode specifications, f and k, for ROUND must be idential with those for the 
calling program; k for ROUND must be one greater than the maximum number of 


digits to be output for any VALUE specified by the calling program. 
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XXKXXXXXXXXKXXKKXXKXKXXKXKXKXKXXKXXKXKX SOURCE PROGRAM XXXXXXXXXXXXXXXKXXXKK KKK KXKXAMK AX 


* 8 8 
€ 1620s FORTRAN IIs COMPUTATION OF WAVE VELOCITIES(FEET/SECOND) AS COMPV O1 
C A FUNCTION OF WATER DEPTH(FEET) AND WAVE PERIOD( SECONDS) e COMPV 02 
E THEORY FROM HeOcoPUBe234(1944)s PROGRAMED BY WeSeWILSON» COMPV 03 
€ JUNE 179 19646 COMPV 04 
DIMENSION CXX (1025) COMPV 05 
TANHF(X) = (EXPF(X)-EXPF(-X)) /(EXPF(X)+EXPF(=X)) COMPV 06 
P = 361415927 COMPV 07 
G = 3202 COMPV 08 
READ 10sNOTT COMPV 09 
10 FORMAT (13) COMPV 10 
DO 1000 NOT=1sNOTT COMPV 11 
READ 209TT COMPV 12 
20 FORMAT (F5el) COMPV 13 
XL = Oe5#G*(TT##200)/(200*P) COMPV 14 
Le xe COMPV 15 
CXXO = TT*G/(2e0*P) COMPV 16 
EEE =| 5ie5 COMPV 17 
BAR = 2e0*P/TT COMPV 18 
DO 2000 K=leL COMPy 19 
DEP = K COMPV 20 
DO 3000 M=1590 COMPyV 21 
CXX(K) = CXXO*TANHF( (BAR¥DEP)/CCC) COMPV 22 
IF (ABSF(CXX(K)-CCC)-e0005) 55300053000 COMPyV 23 
3000 CCC = (CXX(K)+CCC)/200 COMPY 24 
5 IF (SENSE SWITCH 1) 493 COMPy 25 
4 TYPE 900sKsM COMPV 26 
900 FORMAT (2HK=91523H»M=913) COMPyV 27 
3 VALUE = CXX(K) COMPV 28 
CALL ROUND (VALUE9100¢) COMPV 29 
2000 CXX(K) = VALUE COMPV 30 
PUNCH 1009 TTosL COMPV 31 
100 FORMAT (8HPERIOD =sF5e1ls25H SECONDS» MAXIMUM DEPTH =9I597He FEETe)COMPV 32 
PUNCH 200 COMPV 33 
200 FORMAT (/5(5HDEPTHs1X »6HVELCTY 93X)/) COMPV 34 
PUNCH 3009(K»oCXX(K) s9K=1l9L) COMPV 35 
300 FORMAT (5(159F70e293X)) COMPV 36 
PUNCH 700 COMPV 37 
700 FORMAT (///) COMPV 38 
1000 CONTINUE COMPV 39 
TYPE 800 COMPV 40 
800 FORMAT (16HTHIS IS THE ENDe) COMPV 41 
END COMPV 42 

* 8 8 

SUBROUTINE ROUND (VALUE »DEC) ROUND 01 
¢ PROGRAMED BY WeSeWILSON»s JULY 189 19646 ROUND 02 
IVALUE = VALUE * DEC * 106 ROUND 03 
IF (IVALUE) 10051045100 ROUND 04 
104 VALUE = 060 ROUND 05 
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100 


101 


102 
103 


GO TO 103 


JVALUE * VALUE * DEC 


XVALUE = IVALUE 


XVALUE = XVALUE / 1006 


YVALUE = JVALUE 


IF ((XVALUE = YVALUE) = 005) 


VALUE = (YVALUE + le) 


GO TO 103 


VALUE = YVALUE / DEC 


RETURN 
END 


/ DEC 


99.9, 0.0.9.9.0.0.9.9.0.099.9. 0.0.9.2. 0.0.0,0.0.0.0.0.0.00.0.0.04 


1 
400 


XXX XXKXK KKK KK AX XK XK XK XX KK XK KX KKK KXXKXX OUTPUT 


PERIOD = 
DEPTH VELCTY 
1 5060 
6 12283 
ll 16017 
16 18610 
21 19022 
26 19084 
31 20017 
36 20234 


420 SECONDS» 


MAXIMUM DEPTH = 


DEPTH VELCTY 


7082 
13.67 
16064 
18237 
19637 
19293 
20022 
20236 


DEPTH 


10291019101 


INPUT 


ROUND 
ROUND 
ROUND 
ROUND 
ROUND 
ROUND 
ROUND 
ROUND 
ROUND 
ROUND 
ROUND 


KK KH KK HK HK IK KH HK HK MH MM HK IKK HMM HK KK HM KK 


NOTT 


UU 


XM HK MH KKK MA MK KH NHK HH KKK NM HH KK KK KKH KK 


40.0 FEETo 
VELCTY DEPTH 
904§ & 
14040 9 
17.07 14 
18062 19 
19.51 24 
20200 29 
20026 34 
20038 39 
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VELCTY 


10077 
15006 
17045 
18284 
19064 
20007 
20029 
20040 


DEPTH VELCTY 


11287 
15065 
17079 
192004 
19675 
20012 
20032 
20041 


APPENDIX B 
Computer Program for Distribution of Wave Velocity Values 
Over a Grid of Depth Values as a Function of Wave Period 

PROGRAM TITLE: DISTV. 

Variables Used in Program: 

Maiveyetieerien Naver Deriodms (seconds). 

Dodie ceeecee ses NOM a) Sven 4 pOsition) Onvar era GW) aXe lie Vee eaees 
where the grid origin is at X = 0, ¥Y = 0. 

CMAT(I,J)..... For position I,J on the grid, CMAT represents, in the input, 
water depth (feet or fathoms). After conversion, CMAT re- 
presents, in the output, wave velocity (grid units/second). 

MM,NN......... Maximum I and J, respectively, for the grid. 


FMOP.......... Allows conversion of CMAT to feet, if input is in fathoms. 


GRED. se ossee se Grid interval (feet): 
Ingo6 6G OOOO OO 6t 1/2 deep-water wave length (feet) rounded down to the nearest 
integer. 


(CXX(K),K=1,L) For a wave of TT period, CXX is the array of wave velocities 
(feet/second) where K ranges from 1 to L in one-foot incre- 
ments. 

JX,JY...--2--- On each output card JX and JY represents the X and Y grid co- 
ordinates of the first velocity value on that same card. 

Summary Of Program: 

MM, NN, FMOP,GRID,TT, and (CXX(K),K=1,L) are the initial input. Then the 

first row of CMAT depth values is read into the computer. (Fathoms are 

converted if specified by FMOP.) For each depth value the computer "looks 


up" the associated CXX and divides this value by GRID. When all values in 
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the first row have been converted from a depth in feet to a velocity in 
grid units/second, they are rounded to the eighth place to the right of the 
decimal and output. When all rows for a given grid have been similarly 
input, converted, and output, the computer pauses in order that another set 
of data may be input. 

Remarks: 

Each output card includes (at the extreme right-hand side) TT, JX, and JY. 
Depths greater and L are assigned wave velocities equal to those assigned 
for L. The output from this program serves as input for MAIN 1620 and MAIN 
7094 (Appendices D and £), 


MM must be a multiple of 10. 
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XX XK XX KK XXX KK X KKK AXXR XK KAKA KKK AKKX SOURCE PROGRAM XXAXKXKXKKKXAXKAKAAKKKAK AXKAAKAARK AK 


* 8 8 
Cc 16205 FORTRAN IIs DISTRIBUTION OF WAVE VELOCITIES(GRID UNITS/ DISTV 01 
ig SECOND) OVER A GRID OF DEPTHS(FEET OR FATHOMS) AS A FUNCTION DISTV 02 
G OF WAVE PERIOD(SECe) »PROGRAMED BY WeSeWILSON® MAY 19919646 DISTV 03 
DIMENSION CXX(1025) »CMAT (200) DISTV 04 
10 READ 409 MM»oNNoFMOP »GRIDsTT DIstv 05 
40 FORMAT (I149T49F2e09F7e09F501) DISTV 06 
XL = 0¢5*50118%( TT*#*200) DISTV O7 
SG Sede DISTV 08 
READ 209 (CXX(K) sK=loeL) DISTV 09 
20 FORMAT (5(5X9F70293X)) DISTV 10 
PUNCH 100sTT»sGRID DISTv 11 
100 FORMAT (8HPERIOD =9F5e1ls17H SECesGRID SIZE =eF7e0s6H FEETe/) DISTV 12 
DO 3000 J=1sNN DISTV 13 
READ 3059 (CMAT(1I)oI=19MM) DISTV 14 
30 FORMAT (10F4e1) OISTv 15 
IF (FMOP) 29291 DISTV 16 
1 DO 1000 1=19MM DISTV 17 
1000 CMAT(I) = 6e0*CMAT(T) DIsTv 18 
2 DO 2000 I=1 MM DISTv 19 
NVALUE = 1 DISTV 20 
IF (CMAT(1T)) 39200094 DISTV 21 
3 CMAT(T) = -(CMAT(TI)) DISTV 22 
NVALUE = 2 DISTV 23 
4 IF (CMAT(TI)=XL) 69695 DISTV 24 
5 CMAT(I) = XL DISTV 25 
6 XK = CMAT(T) DISTV 26 
K = XK DISTV 27 
CMAT(I) = CXX(K)/GRID DISTV 28 
CALL ROUND (CMAT(I)910¢E8) DISTV 29 
GO TO (200097) sNVALUE DISTV 30 
7 CMAT(I) = =(CMAT(T)) DISTV 31 
2000 CONTINUE DISTV 32 
MM5 = MM/5 DISTv 33 
if al DISTV 34 
BDO 5000 JM=1>5MM5 DISTv 35 
JX = Il-1 DISTV 36 
JY Oo dA DISTV 37 
III = 11+4 DISTV 38 
PUNCH 200s(CMAT(I) sIT2IIs III) sTT»IXsJY DISTV 39 
200 FORMAT (5(F100893X) 92X9FS5elsI40I4) DISTv 40 
Mit ES Shieh teal DISTV 41 
5000 CONTINUE DISTV 42 
3000 CONTINUE : DISTV 43 
TYPE 4009TT DISTV 44 
400 FORMAT (F50e1920H SECeGRID COMPLETEDe) DISTV 45 
PAUSE DISTV 46 
PUNCH 300 DISTV 47 
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300 FORMAT (///) 
GO TO 10 


END 


XX XKK KKK MK KKK KK HM IK KKK KICK XK OOK XK OOKCOK 


20 3 0 

1 5260 

6 120683 

Te LiGrenl 7 

16 18.610 

Aa | WO 

26 19284 

31. 20017 

36 20034 
-00 =-00 -00 
-32 -24 -09 
-00 -00 -00 
-32 -22 008 
-00 -00 -00 
-30 -15 020 


XK KKK KKH KKK HK XXX KKK KK KKK KK XK KKK XKXKKX OUTPUT 


3040. 460 
2 7082 
7 13267 
12 16264 
UY UWISsy 
22 19437 
27 19293 
32 20022 
37 20636 
=-00 -00 -00 -00 
027 030 040 043 
-00 -00 -00 -00 
025 030 039 043 
-00 -00 -00 -00 
030 036 042 045 


-00 
045 
-00 
045 
-00 
049 


PERIOD = 420 SECesGRID SIZE = 3040 
0200000000 0200000000 0200000000 
0200000000 0200000000 0200000000 
-00665132 — 000646053 - 000495395 

© 00671382 000671382 000671382 
0-00000000 0200000000 0200000000 
0200000000 0200000000 0200000000 
- 200665132 -200637171 000473684 

e 00671053 000671382 000671382 
0200000000 0200000000 0200000000 
0.00000000 0200000000 0200000000 
—.00661842 -000585197 000626316 

© 00671382 000671382 000671382 


IN 


9 
14 
17 
18 
19 
20 
20 
20 

41 
055 
43 
056 
-43 
057 


e F 


PUT 

045 4 
040 9 
e207 14 
062 19 
051 24 
200 29 
026 34 
238 39 


EETe 


0200000000 
0200000000 
000655592 
©00671382 
0200000000 
0.00000000 
090649671 
°00671382 
02000009000 
- 200671382 
©00661842 
000671382 
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10677 
15206 
17045 
18284 
19064 
20¢07 
20029 
20¢40 


0200000000 
-000671382 
©00661842 
©00671382 
0200000000 
-000671382 
000661842 
000671382 
0200000000 
- 200671382 
©00669079 
©00671382 


DISTV 48 
DISTV 49 
DISTV 50 


KKK KKM KKK NK KK KK KH KK HK KH KK HI HK IK MK HK KKK 


MM »9NNoFMOP sGRIDsTT 
5 11¢87 
10 15265 
NEY yoy 
20 19204 
AS sos 
30 20¢12 
35 20032 
40 20041 

0000 

1000 

0001 

1001 

0002 

1002 


420 


400 


HK KKK KKK KKM HK HK HK HICH HN IH HHI IK ION IK HK 
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APPENDIX C 
Computer Program to Produce Matrices for Deriving 
Equations of Linear Surfaces 

PROGRAM TITLE: PRMAT. 
Variables Used in Program: 
(X(I),I=1,4).. The X co-ordinates of the four data points. 
(Y(I),I=1,4).. The ¥Y co-ordinates of the four data points. 
SoHMoo600c060d . Matrices used in determiningrthe coefficients for an equation 

of a linear surface of best-fit to data values at four points. 
Summary of Program: 
Let C represent wave velocity at any grid position X,Y. By specifying the 
X,Y values for four positions, this program outputs S and EM. When later 
manipulated with a particular set of C values for these same four positions, 
S and EM produce the coefficients (E's) of linear equation of the form, 
Ca E, + EX + Det This equation represents the least-squares plane of 
best-fit to this particular set of C values. 
Remarks: 
This program was extracted from a general program provided by W. C. Krumbein 
that fits a first, second, or third-order polynomial surface to specified 
regularly-spaced data points by the least-squares method. 
The output from this program serves as spate Oh MAIN 1620 and MAIN 7094 


(Appendices D and £). 
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MX KK KKK KKK KKK KKK XX KK XK KX XXX XMKX SOURCE PROGRAM  XXXKXXXKXKKX KKK KK XK MM HMM KK HHH 


¥*16 6 
Cc 16209 FORTRAN II» PRODUCE MATRICES FOR DERIVING EQUATIONS OF PRMAT O01 
C LINEAR SURFACESe PRMAT 02 
< THEORY FROM WeCeKRUMBEIN» PROGRAMED BY BETTY BENSON» PRMAT 03 
G MODIFIED BY WeSeWILSONs JULY 175 19640 PRMAT 04 
DIMENSION X(4)9 Y(4)9 EM(45s3)9 S(3093) PRMAT 05 
READ 109 (X(I)s9T#194) PRMAT 06 
READ 109 (Y(I)soT=194) PRMAT 07 
10 FORMAT (4F5e0) PRMAT 08 
DO 1000 L=1s4 PRMAT 09 
EM(Lo1) = le PRMAT 10 
EM(L92) = X(L) PRMAT 11 
1000 EM(L»93) = Y(L) PRMAT 12 
DO 2000 [=1s3 PRMAT 13 
DO 2000 J=1s3 PRMAT 14 
S(IsJ) = 000 PRMAT 15 
DO 3000 L=194 PRMAT 16 
3000 S(IsJ) = S(IsJ) +EM(LoT) #EM(L9J) PRMAT 17 
2000 CONTINUE PRMAT 18 
DO 4000 K=193 PRMAT 19 
DIV = S(KsK) PRMAT 20 
S(K9K) = 160 PRMAT 21 
DO 5000 J=153 PRMAT 22 
5000 S(KsJ) = S(KsJ)/DIV PRMAT 23 
DO 4000 I=193 PRMAT 24 
IF (T=K) 19400091 PRMAT 25 
1 DIV = S(ToK) PRMAT 26 
S{I9K) = 000 PRMAT 27 
DO 6000 J=193 PRMAT 28 
6000 S(IsJ) = S(IsJ)I-DIV¥S(K9J) PRMAT 29 
4000 CONTINUE PRMAT 30 
PUNCH 500 PRMAT 31 
500 FORMAT (SOHMATRICES FOR DERIVING EQUATIONS OF LINEAR SURFACES) PRMAT 32 
PUNCH 1005 ((S(TsJ) 9J=193) 912193) PRMAT 33 
100 FORMAT (6F 128) PRMAT 34 
PUNCH 2009 (CEM(LoI) oL=194)sI=193) PRMAT 35 
200 FORMAT (12F6e02) PRMAT 36 
END PRMAT 37 


XH XX KK KK XK KKK KKK KKK MK RK XX KK KXKXKXKKK  TNPUT KXMKKX KKK KKK MH HK KM HN HK KHIM KKM KR KK KK HX 
0.0 1.0 120 0.0 (X( 1) oT=154) 
060 060 160 160 (YCT) oT2=194) 

XX KX XXX KKK KX XXX KK AKK KK KK XK KKKKXKKKKK OUTPUT NXRXXKKKKKKK KK KK MRK KKK HNN NK HK KK HK KKK KK 

MATRICES FOR DERIVING EQUATIONS OF LINEAR SURFACES 

075000000 =e50000000 =-.50000000 =.50000000 1.00000060 0.00000000 


=250000000 0.00000000 14200000000 
1600 1600 1600 100 0200 1-00 1200 0600 0200 0200 1600 1-00 
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APPENDIX D 
Computer Program for Wave Refraction (Linear-Interpolation) 
Using the IBM 1620 

PROGRAM TITLE: MAIN 1620. 

Input Variables: 

S,HM.......... Matrices used in determining coefficients for the plane of 
best-fit to four data points. 

XLABL......... Arbitrary title used to designate each batch of input. 

MM NN. sseece . Maximum values of I and J, respectively, for the CMAT grid 
of vellocity values; T=X +1, J =Y + 1, with the grid origin 
lee@ecech G65 26 = On Mis Oc 

CHECK... .ccces Allows either a two-dimensional CMAT to be input, or a single- 
rowed CMAT (extending from shore to deep water) to be input 
and all other rows to be made identical to the first. (This 
procedure creates a set of grid values which can be character- 
ized by straight and parallel contours. ) 

TBs Sooo .sss.. Wave period (seconds). 

NOJ.......---- Number of wave rays to be run in each batch. 

D..eeeeeceeese- Distance incremented between successive points along a ray 
path. 

CMAT(I,J)..... Wave velocity (grid fiateyeccen a) at grid position I,J. 

A..ecccseeeees Angle (degrees) measured from the direction of increasing X 
along the X-axis moving counter-clockwise to the dire¢tion of 
travel of a wave ray. 


X,Y........-.- Grid origin position for a wave ray. 
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Output Variables: 


WIWABIS Goo G0000 Defined previously. 

TR h9 OltiOIOO.G Defined previously. 

INOIWs 6660 see... Davch number. 

Woosodana00000 Ray number for batch NOT. 

MIX cocaod0000 - Number assigned to a point along a ray path where calculations 


are made. MAX = 1 at the origin point of ray. 

FOOL, WAR 506000 X,Y co-ordinates for a given MAX. 

ANGLE......... Angle (degrees) between X-axis and ray, as previously defined 
for A. 

PCTDIF.....+-. Percent difference between the value of one of the four 
points to which a plane is fit and that corresponding value 
of the plane (see text, "Interpolation Surfaces"). 

Variables in Common: 

S, HM.......-. Defined previously. 

Heesceeeeeeees Coefficients of the equation of a plane of best-fit to four 
grid points. 

YVW......-.... Matrix used in determining E's. 

CWE gc00G0006 Defined previously. 

Crccccceeeeeee Values of the four CMAT values to which a plane is fit. 

XLABL......... Defined previously. 


D.....--2----- Defined previously. 


IM6 60000000 -.-. Defined previously. 
QiMooac0o00000 Wave-velocity for a given MAX. 
Io aga00000000 Number of iterations used in obtaining the curvature of a 


wave ray for a given MAX. 


NGO........... Designates whether a wave ray has moved within one grid unit 
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of the edge of the grid. 

AMM, ANN....... Used in determining NGO. 

MAX......-.-.. Defined previously. 

MUEPG 6 od0000 -.. Designates whether the last two curvature estimates for a 
given MAX are less than 0.00009/D, whether the 18th and the 
20th estimates are less than this value, or whether neither 
of the above is true. 

Summary of Program: 

MAIN reads S and EM and sets NOT = 1 to denote batch one. It then reads the 

information (XLABL, MM, NN, CHECK, TT, NOT, D, and CMAT) for processing the 

first batch of rays. For N = 1 the program reads the information (Gi, 35° ekavel 

A) for processing the first wave ray. MAIN converts A from degrees to radians 

and then punches the title information (XLABL, TT, NOT, N, and the column 

headings). Control is then transferred to the RAYN subroutine. When RAYN 
has determined the path of the first ray, control is transferred back to MAIN. 

MAIN then reads the information for the second ray (N = 2) and proceeds as 

before. When NOJ ray paths have been determined for NOT = 1, the computer 

pauses allowing time to load the necessary information for the second batch. 

After pressing start, MAIN sets NOT = 2, and the information for this batch 

is read. MAIN then continues as before. The program is terminated when the 

desired number of batches has been processed. 

Remarks: 

5; 
This program calls four subroutine/RAYN, SURFCE, MOVE, and ROUND. Descrip- 


tions of these subroutines follow. (ROUND has been described previously. ) 
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SUBROUTINE TITLE: RAYN. 

Summary of Subroutine: 

RAYN calls SURFCE to obtain FK (ray curvature) for MAX = 1. MAX is then set 
equal to two. RAYN calls MOVE in order that X, Y, A, and FK be obtained for 
MAX = 2. ANGLE (equal to A but expressed in degrees instead of radians) and 
PCTDIF are then computed. XXX and YYY are set equal to X and Y, respectively, 
so that significance is not lost when the values to be output are rounded. 
ROUND is called for XXX, YYY, ANGLE, and PCTDIF before they are punched. MAX 
is incremented by one, and the entire process continues until 1) MAX = 500 , 
2). CXY is' equal to or less than zero, 3) MIT = 3, or 4) NGO = 3. The first 
condition guards against a ray going in an endless circle. The second pre- 
vents a ray from leaving the water and going over land. The third prevents 
invalid curvature approximations. The fourth condition stops a ray if it 
reaches the edge of the grid. Any of these four conditions causes the 
message, "RAY STOPPED" to be punched and transfers control back to MAIN. 

If MIT = 2, the message "CURVATURE APPROXIMATED" is punched and computations 
proceed as usual. 

Remarks: 

The use of Sense Switch 3 keeps the operator informed of the current value 
of MAX. 

SUBROUTINE TITLE: SURFCE. 

Variables Used in Subroutine: 

I,FI......+.-- X + 1 rounded down to the nearest integer. 

J,FJe.e+eee--- Y + 1 rounded down to the nearest integer. 

Wbooaocco00000 o& we dh o Wilko 


Wabooccaccc0000 er do Wo 
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We UN 00000000 Values of FI and FJ, respectively, the last time SURFCE was 
called. (Not valid for MAX = 1.) 

IMK&5 g00000G0000 Curvature of the wave ray at X,Y. 

Summary of Subroutine: 

For a specified A and X,Y position on the CMAT grid, SURFCE first determines 

I, J, XL, and YL. If MAX does not equal one, SURFCE checks to see if ZI = 

FI and ZJ =FJ. If both equalities are true, the # values from the previous 

operation of SURFCE are still valid, and the CXY and FK computations can be 

made directly. Otherwise, it is necessary to derive new E values, using C, 

IM, and S, before computing CXY and FK. 

Remarks: 

Program steps corresponding to SURFCE2O0 through SURCH2/7 in the listing that 

follows were obtained from a program provided by Dr. W. C. Krumbein. (This 

program is described with the PRMAT program in Appendix Gs) 

SUBROUTINE TITLE: MOVE. 

Variables Used in Subroutine: 

FKBAR......-.-- Curvature used to obtain DELA for a given IT. 


TMM 6 obcGoad MINS yor ITU = AUS). 


ll 


FKKP.......... FKBAR for IT = n-l where current FKBAR is for IT =n. (For 
IT = 1, FKKP equals the FKBAR used in determining XX and YY 


the previous time MOVE was called.) 
PKK .cc-ceee-- X, Y, AS and HK wailues, respectively.) for) Max snr iewinercce 


MAX = n when MOVE was called. 


DELX. 2. ccee ~. XX = X. 
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Summary of Subroutine: 

If MAX equals three or more when MOVE is called, FKBAR for the previous MAX 
is used in obtaining approximations of XX, YY, and AA. (For MAX = 2, FKBAR 
is set equal to FK.) SURFCE is called and returns FKK for this approximation 
at XX,YY. KBAR is then redefined as (FK+FKK)/2. If the difference between 
FKBAR and FKKP is less than 0..00009/D, the current XX, YY, AA, and FKK values 
are accepted for the new point. If the difference is greater, FKKP is set 
equal to FKBAR, and the current FKBAR is used to obtain another set of XX, 
YY, AA, and FKK values. The difference between FKBAR and FKKP is again 
tested. This cycle may repeat a maximum of 20 times before termination. If 
the cycle stops before IT = 20, MET is set equal to one. If the cycle stops: 
at IT = 20, and if the difference between FKBAR and FKKPP is less than 
0.00009/D, then MIT is set equal to two, and FKBAR is defined as (FKBAR+FKKP)/2 
for obtaining XX, YY, AAA, and FKK. If IT = 20, and this difference is 
greater than 0.00009/D, then MIT is set equal to three. When MIT = 3, 
control is transferred back to RAYN immediately. When MIT = 1 or 2, XX and 
YY are tested to see if the new point has reached the édge of the grid. 


NGO = 2 if the ray has reached the edge of the grid, and NGO = 1 if it has 


not. 

Remarks: 

MIT = 1 when the curvature approximations are converging to a single vaiiue. 
MIT = 2 when the approximations are oscillating between two values. In this 


case the average of the two values is taken as the new curvature. MIT = 3 
when the approximations are oscillating among three or more values. No 

valid curvature approximation can be made in this case; this is one basis for 
the terminiation of a ray. 


The use of Sense Switch 2 allows the operator to observe successive IT and 
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FKBAR values. 
The maximum difference (0.00009/D) in the curvature approximations is such 


that the output ANGLE is significant in the hundredths digit. 
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KX KH KKK MK KKM HHH HMM HK KK HK HK KKK HK HK KK 


NA * 


403 FORMAT (///12A4/8HPERIOD =9F50e1s6H SECes910H BATCH NOceI399H» RAY 
1391He//4X 9 3HMAX 9 6X 9 LHX 9 8X 9 LHY 9 8X9 SHANGLE 9 4X 9 6GHPCTDIF 9 
2//1792F9e29Fl1lle2) 


15 


3 


INOec» 


1620s FORTRAN IIs LINEAR-INTERPOLATION WAVE REFRACTION PROGRAMe 

ORIGINAL PROGRAM BY GRISWOLDsNAGLEsAND MEHRe THIS PROGRAM ADAPTED 
FROM ORIGINAL BY WILSON»HARRISON »sKRUMBEINs AND BENSONe 
DIMENSION S(393) 5EM( 493) 5E(3) »¥YVW(3) sCMAT (50951) 9C(4) »XLABL(12) 
COMMON SsEMsEsYVWsCMAT sC oXLABL »DoTT »CXY 9IT 9NGO 9 AMMo ANN 9 MAX o MIT 


READ 59(¢€S(I9J) »J=193) s1=193) 
FORMAT (6F1208) 

READ 7o((CEM(LoI) oL=1904) »I=193) 
FORMAT (12F 6602) 

NOT = 1 

READ 4009 XLABL 

FORMAT( 12A4) 

READ 402 9MMoNN» CHECK 9 TT oNOJ 9D 
FORMAT (2149F3ce097X9F50e19159F40l) 
AMM = MM=1 

ANN = NNe1l 

IF (CHECK) 10910920 

READ 119((CMAT(I 9J) »IT=19MM) »J=1 9NN) 
FORMAT (5(F100893X)) 

GO TO 14 

Jzil 

READ lls (CMAT(I09J) 9I=159MM) 

DO 77 J=2o9NN 

DO 77 I=1l9MM 

CMAT(T9J) = CMAT(T91) 

DO 15 N=l»5 NOJ 

READ 69A0XsY 

FORMAT (F70e292F6e02) 

MAX = 1 

PUNCH 4039XLABL»oTT  »NOT oN oMAXoX9YOA 


A=A¥ 00174532925 
CALL RAYN (X9YsA) 
CONTINUE 

PAUSE 

NOT = NOT + 1 

GO TO 9998 

END 


SUBROUTINE RAYN (XsY9A) 


DIMENSION S(393) 9EM( 493) 9E(3) sYVW(3) »CMAT(50951) 9C(4) » XLABL(12) 
COMMON SoEMsEsYVWsCMAT 9C »XLABL sDoTT sCXYsIT #»NGO» AMM 9 ANN 9 MAX o MIT 


CALL SURFCE (XsYsAsFK) 
MAX=1+MAX 
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MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 
MAIN 


RAYN 
RAYN 
RAYN 
RAYN 
RAY¥N 
RAYN 


SOURCE PROGRAM  XXXXXKKXXXKXXKKMAKK MAKKAH KM KK KX HA KK NK 


100 
103 
101 
399 


396 
395 
200 
397 


12 


13 


8 6 


MPN 


318 


IF (SENSE SWITCH 3) 1002101 

TYPE 103 9MAX 

FORMAT (4HMAX=914) 

IF (MAX=500) 399915915 

CALL MOVE (XsYoAsFK) 

TF -CEXY) 2159159396 

GO TO (3979395515)s MIT 

WRITE 2009 MAX 

FORMAT (32HCURVATURE APPROXIMATED FOR MAX =914) 
ANGLE=A*57e29577951 

XXX = X 

WOAL Evg 

PGTDINE = ABSE CUG(3) -E CRE C2). =F 3) G03) 21 0 Oe 
CALL ROUND (XXX9100e) 

CALL ROUND (YYY:100e) 

CALL ROUND (ANGLE »1000e) 

CALL ROUND (PCTDIF 9100) 

PUNCH 129MAX9XXX9YYY o9ANGLE sPCTDIF 
FORMAT (I1792F9e29Flle29F10e0e1) 

GO TO (3915) »NGO 

PUNCH 13 

FORMAT (12HRAY STOPPEDe) 

RETURN 

END 


SUBROUTINE SURFCE (XsYsAo9FK) 


DIMENSION S(393) 2EM(493) E(3) sYVW(3) sCMAT (50951) 9C(4) sXLABL(12) 
COMMON SoEMoEsYVWsCMAT »CsXLABL »DoTT »CXY9IT»NGOsAMM oANN 9 MAX o MIT 
I1=X+lo 

J=Y+1le 

FI=I 

FJ=J 

XL=X+1le-FI 

YL=Y+1le-FJ 

IF (MAX=1) 19194 

VE CZT—-F1) Lsi2is1 

WE (ZI=F)) 15391 

Zao 

Zj) SB Fd) 

C(1)=CMAT(I9J) 

C(2)=CMAT(I4+1sJ) 
C(3)=CMAT(I4+19J+1) 
C(4)=CMAT( I 9J+1) 

DO 318 Il=193 

YVW(TT) = Oo 

DO 318 L=1s4 

YVW(IT) = YVWOITT)4C(L)*EM( Lo IT) 
DO 319 II=193 

ECII) = Oc 
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RAYN O7 
RAYN 08 
RAYN 09 
RAYN 10 
RAYN 11 
RAYN 12 
RAYN 13 
RAYN 14 
RAYN 15 
RAYN 16 
RAYN 17 
RAYN 18 
RAYN 19 
RAYN 20 
RAYN 21 
RAYN 22 
RAYN 23 
RAYN 24 
RAYN 25 
RAYN 26 
RAYN 27 
RAYN 28 
RAYN 29 
RAYN 30 


SURFCEO1 
SURFCE02 
SURF CE03 
SURFCE04 
SURFCE05 
SURFCE06 
SURFCEO7 
SURFCEO8 
SURFCEO9 
SURFCE10 
SURFCE11 
SURFCE12 
SURFCE13 
SURFCE14 
SURFCE15 
SURFCE16 
SURFCE17 
SURFCE18 
SURFCE19 
SURFCE20 
SURFCE21 
SURFCE22 
SURFCE23 
SURFCE24 
SURFCE25 


319 


102 
104 
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On Ww 
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DO 319 JJ=193 
ECII) = ECIIV+S(TIsJJ)*#YVW( JJ) 


CXY = E(1) + E(2)¥XL + EU3)*YL 

FK =(E(2)*SINF(A)-E(3)*COSF(A)I/CXY 
RETURN 

END 


SUBROUTINE MOVE (XsYsAo0FK) 


DIMENSION S(393) 0EM( 4093) 9E(3) sYVW(3) sCMAT(50951) 0€{4) oXLABL(12) 
COMMON SoEMoEsYVWoCMAT oC sXLABL sDs TT sCXY oI TsNGOoAMMoANR oMAX oMIT 
IF (MAX = 2) 10291029104 

FKBAR=FK 

MIT = 1 

DO 20 IT#1920 

DELA=FKBAR*D 

AA=A+DELA 

ABAR=A+e5*DELA 

DELX=D#COSF (ABAR) 
DELY=D#SINF(ABAR) 

XX=X+DELX 

YY=Y+DELY 

GO TO (10196) MIT 

CALL SURFCE (XX9YYs2AAs9FKK) 

IF (CXY) 38938910 

FKBAR = 00e5 *® (FK + FKK) 

IF (SENSE SWITCH 2) 8989899 

TYPE 9009 ITsFKBAR 

FORMAT (149E14-8) 

IF (IT = 18) 593799 

FKKPP = FKBAR 

TF (MAX = 2) 79709 

IF (IT = 1) 2092009 

IF (ABSF(FKKP=FKBAR) = (0-00009/D)) 696920 
FKKP = FKBAR 

IF (ABSF(FKKPP = FKBAR) = (0200009/0)) 18518917 
MIT = 3 

GO TO 38 

FKBAR = 0065 * (FKBAR + FKKP) 

MIT = 2 

GO TO 39 

NGO = 1 

IF ((XX—100)*( (AMM—10e00)—XX))29293 
IF (€YY—1¢0)*( (ANN=100)-YY) 129298 
NGO = 2 

X = XxX 

Y 2 YY 

A = AA 

FK = FKK 

RETURN 

END 
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SURFCE26 


SURFCE27 
SURFCE28 
SURFCE29 
SURFCE30 
SURFCE31 
MOVE 01 
MOVE 02 
MOVE 03 
MOVE 04 
MOVE 05 
MOVE 06 
MOVE 07 
MOVE 08 
MOVE 09 
MOVE 10 
MOVE 11 
MOVE 12 
MOVE 13 
MOVE 14 
MOVE 15 
MOVE 16 
MOVE 17 
MOVE 18 
MOVE 19 
MOVE 20 
MOVE 21 
MOVE 22 
MOVE 23 
MOVE 24 
MOVE 25 
MOVE 26 
MOVE 27 
MOVE 28 
MOVE 29 
MOVE 30 
MOVE 31 
MOVE 32 
MOVE 33 
MOVE 34 
MOVE 35 
MOVE 36 
MOVE 37 
MOVE 38 
MOVE 39 
MOVE 40 
MOVE 41 
MOVE 42 
MOVE 43 
MOVE 44 


XXXXXXXXXXXXXXXXXXXK KX KK KK XX KK KK KKK 


275000000 =.50000000 
=250000000 000000000 
NoGO WMEOO AGOO ” W5OO 

1620s GRID OFF VAeCAPES» 
Ao) oe © 400 
0200000000 0200000000 
0200000000 0200000000 
— 200665132 - 200646053 
© 00671382 000671382 
0200000000 0200000000 
0200000000 0200000000 
—.00665132 -200637171 
©00671053 000671382 
0-00000000 0200000000 
0.00000000 0.00000000 
- 00661842 -.00585197 
2° 00671382 000671382 
0200000000 0200000000 
0200000000 0200090000 
= 000646053 000354276 
200671382 000671382 
0200000000 0200000000 
0200000000 0200000000 
-2 00612500 ©00514803 
© 00671382 ©00671382 
0200000000 0200000000 
0200000000 0200000000 
= 000495395 000637171 
200671382 200671382 
0200000000 0200000000 
0200000000 -200671382 
200422039 ©00660197 
©00671382 ©00671382 
0200000000 0200000000 
0200000000 -200669737 
© 00531908 200666447 
©00671382 ©00671382 
0200000000 0200000000 
0200000000 - 200669737 
© 00655592 000666447 
© 00671382 200671382 
000000000 0200000000 
-.00669737 000666447 
200663487 000668421 
200671382 000671382 
0200000000 0200000000 
=. 00669079 —000666447 
© 00665132 000668421 
© 00668421 000671382 


INPUT 


-~.50000000 -.50000000 
1200000000 
0.00 1200 1-00 0-00 
LINEAR INTERPos 8/5/646 
3 05 
0e00000000 0-00000000 
0600000000 0.00000000 
= 200495395 200655592 
200671382 200671382 
0e00000000 0600000000 
0200000000 0.00000000 
000473684 2©00649671 
200671382 200671382 
0600000000 0.00000000 
0200000000 -.00671382 
© 00626316 200661842 
200671382 200671382 
9200000000 0.00000000 
-.00671382 -.00671382 
000657895 200669079 
200671382 2©00671382 
0e00000000 0.00000000 
-000671382 -e00670395 
200665132 200669737 
000671382 200671382 
0e00000000 0200000000 
=200671382 -—.00670395 
2000661842 200671382 
200671382 200671382 
0200000000 0.00000000 
—000670395 200666447 
200669079 °00671382 
200671382 °00671382 
02000000000 0.00000000 
=000669737 —000666447 
200670395 200671382 
200671382 200671382 
0200000000 0.00000000 
~.00668421 -.e00665132 
200669737 200671382 
200671382 200671382 
0e00000000 0200000900 
-000665132 000652632 
200669737 200671382 
200671382 200671382 
0e00000000 0.00090000 
=000665132 -000619737 
2000670395 200671382 
000671382 200671382 
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0200 


1200000000 


0200 


0e00000000 
-°00671382 
©00661842 
°©00671382 
0200000000 
-—2©00671382 
©00661842 
200671382 
0e00000000 
- 200671382 
°00669079 
©00671382 
0200000000 
-— 200669079 
000671382 
©00671382 
0200000000 
— 000666447 
©00671382 
©00671382 
0200000000 
—200661842 
©00671382 
000671382 
000000000 
-200649671 
¢00671382 
©00671382 
0200000000 
-—e00612500 
©00671382 
©00671382 
0200000000 
— 000495395 
°00671382 
000671382 
0200000000 
°00390461 
©00671382 
©00671382 
0200000000 
©00531908 
©00671382 
©00671382 


0.00000000 


1-00 1-600 


40 


FPP ELE 
eee 8 @ © e@ Co 


9o00O0 00000 


> 
e 


XXXKKXKXKKXKK RX KKK KX XX KX KX KKK KKK KAKA KKK KK 


Ll 


0-00000000 
— 200667434 
© 00665132 
e 00671382 
0200000000 
- 000666447 
© 00665132 
© 00669737 
0200000000 
— 200665132 
© 00667434 
e 00669737 
0200000000 
- 000655592 
© 00665132 
e00669079 
0-00000000 
-—2°00649671 
e00665132 
°00671382 
0-00000000 
— 200646053 
© 00661842 
© 00671382 
0-00000000 
-—-00649671 
© 00665132 
e 00671382 
0.00000000 
-—200649671 
© 00667434 
©00671382 
0200000000 
-—¢00649671 
° 00665132 
© 00671382 
0200000000 
-— 200641776 
e 00663487 
© 00671382 
0200000000 
-200637171 
e 00665132 
© 00671382 


12060 14650 
120060 15036 
120¢0 16622 


0200000000 0200000000 0200000000 0200000000 4.0 0 11 
— 200666447 -000652632 -000547368 ©00649671 400 Bait 
000666447 000669737 ©00671382 ©00669079 4-0 10 11 
©00671382 ©00671382 e00671382 e00671382 400 15 11 
0200000000 0e00000000 000000000 - 00669079 400 ahr 
- 000665132 -000649671 ©00257237 ©00657895 400 Sel 
000667434 ©00671053 ©00671382 ©00669079 400 10 12 
©00671382 ©00671382 ©00671382 ©00671382 GeO) eo 2 
0200000000 0e00000000 0200000000 —000666447 40 @). - 342} 
-200661842 -— 000626316 ©00514803 ©00661842 400 Gy AN) 
000666447 ©00669079 © 00671382 ©00669737 4e0 10 13 
000671053 ©00671382 © 00671382 ©00671382 4¢0 15 13 
0200000000 02e00000000 0200000000 — 000669737 40 Oo 14 
-—000641776 -.00585197 ©00637171 000661842 4.0 5 14 
000667434 ©00670395 000671382 000671382 4e¢0 10 14 
000671053 ©00671382 © 00671382 ©00671382 4¢0 15 14 
0200000000 0200000000 -—200668421 -000661842 420 @) 3G) 
-200637171 ~-2©00257237 e 00646053 ©00657895 400 5S 
000667434 °00669079 «00671053 ©00671382 460 10 15 
©00671382 000671382 200671382 ©00671382 40.0 15 15 
0200000000 0200000000 - 200665132 -—000657895 4.0 Oo 16 
- 000626316 000514803 © 00655592 °00657895 400 3 Ue 
000667434 ©00669079 © 00671382 ©00671382 4e-0 10 16 
©00671382 e00671382 ©00671382 °00671382 4-0 15 16 
0200000000 0200000000 -— 200663487 —000657895 420 Oo 17 
— 000626316 ©00547368 °©00660197 °00660197 400 Sli 
\ 200669079 ©00670395 °00671382 ©00671382 4e0, 10 17 
000671382 00671382 °00671382 °00671382 4e0\ 15 17 
0200000000 000000000 -— 200660197 -—200657895 4.0 @) alts) 
-\.00585197 ©00637171 © 00652632 ©00661842 \4e0 Sees 
000668421 000669737 ©00670395 ©00671382 400 10 18 
600671382 ©00671382 e00671382 ©00671382 ‘460 \15 18 
0200000000 0200000000 -200657895 =000655592 4e0 @) oO) 
= 200422039 ©00649671 °00652632 °00660197 420 By 1S) 
000666447 ©00671053 ©00671382 °00671382 42.0 \10 19 
000671382 000671382 °00671382 ©00671382 400 15 | 19 
0200000000 — 000663487 - 200655592 = 000652632 4-0 0 \20 
0200000000 000655592 ©00655592 ©00655592 40 5 \20 
e00665132 °00669079 °00671382 ©00671382 400 10 (20 
©00671382 000671382 ©00671382 °00671382 400 15 20 
0200000000 -—200660197 —200655592 200649671 420 Oo 21 
000547368 000655592 °00655592 °00655592 40 5) 2.1! 
000667434 ©00670395 ©00671382 000671382 400 10 21 
000671382 °00671382 ©00671382 ©00671382 4e0 15 21 
0250 RAY 1 
02099 RAY 2 
03049 RAY 3 


S| 


XXXXXXXXXXKXXXXXKXXXXXXXKXKXXXXXXKXKXXXKXKX OUTPUT XXXXXXXXXXXKXKXKXKXKXKXKXKXKXKKK KKK KKK KKK KK 


16205 GRID OFF VAeCAPESs LINEAR INTERPe.98/16/64e 


PERIOD = 4¢0 SECes BATCH NOw 1s RAY NOe le 
MAX X Y ANGLE PCTNIF 
] 14650 2050 120-00 
2 14025 2093 120207 el 
3} 14200 3037 120614 00 
4 13275 3080 120022 020 
5 13250 4023 12028 el 
6 13024 4066 120.23 el 
U 12.99 5099 120¢51 3 
8 12¢74 Dei2 12080 °3 
9 1248 5295 121.210 03 
10 12022 6028 121228 0-0 
11 11-96 608) 121048 02 
2 11270 Te23 121.270 020 
13 11643 7266 121283 0.0 
14 Wiloahy 808 121294 ol 
WS) 19290 8251 122019 02 
16 10.64 8293 122059 e2 
Wy NOSS7 9035 12288 el 
18 19.09 GeTl 123295 el 
19 9281 10.18 125027 404 
20 9051 10.58 129256 44 
21 9-18 10095 133285 4o4 
22 8.75 11222 162047 3023 


RAY STOPPED. 


1620s GRID OFF VAeCAPESs LINEAR INTERPe 98/16/64e 


PERIOD = 4e0 SECes BATCH NOe ls RAY NOco 26 

MAX x Y ANGLE PCTDIF 
1 1536 2099 120200 
2 MB) aba 3042 120-00 000 
3 146286 3086 120-00 0.0 
4 14.61 4079 120-00 0-0 
5 14036 4072 120200 000 
6 14.11 52016 120200 000 
U 13286 5059 120-90 0e0 
8 13661 6092 12000 020 
9 136026 6045 120-00 0560 — 
10 13611 6089 120200 00 
11 12286 7032 120.02 000 
2. 1261 T0775 120-06 0-0 
13 1236 8019 120e11 020 
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14 12611 8e62 120617 020 
15 11.286 9205 120224 020 
16 11.60 9248 120230 0e0 
W7/ 1135 9.91 120237 0e0 
18 11.10 1034 120044 000 
19 10.85 1078 12051 el 
20 10.59 11.21 120057 000 
2 10034 110264 12065 020 
22 10208 12207 120671 el 
23 9.83 12250 120.88 ol 
24 9057 1292 121019 el 
25 9e31 13035 121641 el 
26 9205 12278 121055 ol 
27 8.78 1420 121-98 05 
28 8052 14063 122672 05 
29 8024 15604 123027 04 
30 7091 15042 139233 2961 
3)il 7046 15263 W/o als} 29] 


RAY STOPPEDe 


1620s GRID OFF VAeCAPESs LINEAR INTERPe 98/16/64e 


PERIOD = 4ce0 SECes BATCH NOs 1s RAY NOc 36 
MAX x Y ANGLE PCTDIF 
1 16022 3049 120.200 
Q 15297 2092 120200 000 
3 15072 4236 120200 020 
4 15047 419 120200 0-0 
5 15022 5022 120.00 020 
6 14-97 5266 120200 000 
v 14.72 6009 120-00 000 
8 14047 6052 120200 000 
9 14622 6095 120200 0-0 
10 13097 72039 120200 0-0 
i 13672 7082 120-00 020 
UZ 13047 8225 120200 0-0 
LES 13622 8-69 120-00 000 
14 12697 9012 120-03 00 
15 12272 9255 12008 000 
16 12047 9298 120014 0-0 
17 12022 10642 120019 020 
18 11.97 10285 120024 00 
19 1le71l 1128 12035 020 
20 11046 11.71 120651 000 
21 11621 12014 12062 000 
22 10095 12057 120.68 el 
23 10.70 13-00 120671 el 
24 10044 13243 12072 ol 
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25) 10019 13086 
26 9293 1429 
2at 9067 14672 
28 9e41 15015 
29 9216 15258 
30 8290 16200 
Bil 8.63 16043 
32 837 16086 
3)3) Bell 1728 
34 7284 17270 
35 7056 18.11 
36 7026 18.252 
37 6091 18.87 


RAY STOPPED. 


12073 
120-81 
120-96 
121212 
12129 
12144 
121-56 
121268 
121.79 
123651 
125042 
125092 
144.59 
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APPENDIX E 
Computer Program for Wave Refraction (Linear-Interpolation) 
Using the IBM 7094 
PROGRAM TITLE: MAIN 7094. 
Summary of Program: 
This program is essentially the same as that of MAIN 1620; both programs 
give identical results if processing the same input data. The chief differ- 
ences between the programs are as follows: 
(1) MAIN 7094 is written in FORTRAN IV, while MAIN 1620 is written in 
FORTRAN IT. 
(2) NOTT (the maximum number of batches to be run) is specified for MAIN 
7091. 
(3) Sense Switches are not used in MAIN 709}. 
(4) ROUND subroutine is not used in MAIN 7094 since the 7094 computer auto- 
matically rounds off output data. 
(5) Due to the large storage capacity of the 7094, dimensions for CMAT in 
MAIN 7094 may greatly exceed those in MAIN 1620. 
(6) Since the 7094 computer outputs to a printer, MAIN 7094 has been pro- 
grammed to print title information at the top of each output page. 
(7) MAIN 7094 operates approximately 200 times faster than MAIN 1620. 
Remarks: 
MAIN 1620 is best suited for experimental work connected with further pro- 
gram refinements and development, while MAIN 7094 is best suited for pro- 


cessing large numbers of rays. 
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XXXMXKXKXKXKX KKK KX XXKXKKKXKXKXKKXXKXXX SOURCE PROGRAM  XXXXXKXXKXKK KKK KM KKK HK KKK KKK KKK KK 


$IBFTC MAIN MAIN O1 

G 70945 FORTRAN IVs LINEAR-INTERPOLATION WAVE REFRACTION PROGRAMe MAIN 02 

C ORIGINAL PROGRAM BY GRISWOLD»sNAGLE»AND MEHRe THIS PROGRAM ADAPTED MAIN 03 

Cc FROM ORIGINAL BY WILSON»HARRISON»KRUMBEINS AND BENSONe 8/4/640 MAIN 04 

DIMENSION S(393)9EM( 493) 9E(3) sYVW(3) »CMAT(100982) 9C(4) »XLABL(12) MAIN 05 

COMMON 5 » EM o — » YVW » CMAT iG MAIN 06 

COMMON XLABL » D 5 Tt » CXY OT » NGO MAIN O07 

COMMON AMM » ANN » MAX » NOT » N »9 MIT MAIN 08 

READ (595) ((S(Is9J)9J=193)s9l=193) MAIN 09 

5 FORMAT(6F1268) MAIN 10 

READ (597) ((EM(LsI)9L=194) sl=193) MAIN 11 

7 FORMAT(12F602) MAIN 12 

READ (59500) NOTT MAIN 13 

500 FORMAT (13) MAIN 14 

DO 399 NOT=1leNOTT MAIN 15 

READ (59400) XLABL MAIN 16 

400 FORMAT (12A6) MAIN 17 

READ (59402) MMsNN»CHECK »TTsNOJsD MAIN 18 

402 FORMAT (2149F30007XsF5elsI59F4ol) MAIN 19 

AMM = MM=1 MAIN 20 

ANN = NN-1 MAIN 21 

ITF (CHECK) 10910920 MAIN 22 

10 READ (5911) ((CMAT(I9J)9T=19MM) 9J=19NN) MAIN 23 
11 FORMAT (5(F100¢8s3X)) MAIN 24 | 

GO TO 14 MAIN 25 

AG) Je i MAIN 26 

READ (5911) (CMAT(19J)9I=19MM) MAIN 27 

DO 77 J=29NN MAIN 28 

DO 77 I=19MM MAIN 29 

77 CMAT(T9J) = CMAT(I91) MAIN 30 

14 DO 15 N=l»s NOJ MAIN 31 

READ (596) AsXsY MAIN 32 

6 FORMAT (FT70e292F6e02) MAIN 33 

MAX = 1 MAIN 34 

WRITE (69403) XLABL»oTTsNOToNoMAX 9X 9A MAIN 35 

403 FORMAT (1H191246/9H PERIOD =9F5016H SECess1l0H BATCH NOcosl399H» RAMAIN 36 

LY NOes1391He//4X 9 3HMAX 9 6X 9 LHX 98X 9 LHY 9 8X 9 SHANGLE 94X 96HPCTDIF// MAIN 37 

21792F9e29Flle2) MAIN 38 

A=A% 00174532925 MAIN 39 

CALL RAYN (X9YsA) MAIN 40 

15 CONTINUE MAIN 41 

399 CONTINUE MAIN 42 

WRITE (699999) MAIN 43 

9999 FORMAT (18H1 THIS IS THE ENDe) MAIN 44 

CALL EXIT MAIN 45 

STOP MAIN 46 

END MAIN 47 
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C 


SIBFTC RAYN 


G 


SUBROUTINE RAYN (X9YosA) 


DIMENSION S(393)9EM( 493) 9E(3) sYVW(3) sCMAT( 100982) 9C(4) sXLABL(12) 


COMMON S 9 EM 9 |e » YVW 9 CMAT 
COMMON XLABL » D » TT » CXY OH) 0 
COMMON AMM » ANN » MAX » NOT » N 
CALL SURFCE (X9YsAsFK) 

MAX=1+MAX 


IF (MAX = 800) 399915915 

CALL MOVE (XeYsAsFK) 

IF (CXY) 159159396 

GO. TO (3979395915)s5 MIT 

PUNCH 2005 MAX 

FORMAT(33H CURVATURE APPROXIMATED FOR MAX =914) 
ANGLE=A*57029577951 

IF (MOD(MAXs40)) 2025920 

WRITE (697) XLABLsTT»NOT9N 


FORMAT (1H1912A6/9H PERIOD =9F5e1l96H SECess10H BATCH NOes1399H»s 
1Y NOesI391He//4X 9 3HMAX 9 6X 9 LHX 9 8X 9 LHY 98X 9 SHANGLE 9 4X s6HPCTDIF//) 


PCTDIF = ABSF((C(3)-E(1)-E(2)-E(3))/C(3))*1006 
WRITE (6912) MAX sX»9Y ANGLE sPCTDIF 

FORMAT (1792F9e2sFlle2sFl10el) 

GO TO (3915) »NGO 

WRITE (6913) 

FORMAT (13H RAY STOPPEDe) 

RETURN 

END 


S$IBFTC SURFCE 


PN 


SUBROUTINE SURFCE (XsYsAsFK) 


DIMENSION S$(333)59EM( 493) 9E (3) sYVW(3) »CMAT(100982) 9C(4) »XLABL(12) 


COMMON §S » EM 9 E » YVW » CMAT 
COMMON XLABL » D 9 TT » CXY » IT 
COMMON AMM » ANN » MAX » NOT » N 
I=X+1l. 

J=Y+1le 

FI=I 

FJ=J 

XL=X+1le-FI 

YL=Y+1le-FJ 


IF (MAX=-1) lslo4 


Ife (Zur) oeow 
IF (ZJ-FJ) 19391 
Zu S [py 
ZJ = FJ 


C(1)=CMAT(T J) 
C€(2)=CMAT(I+19J) 
C(3)=CMAT(I+19J+1) 
C(4)=CMAT( I 9J+1) 
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9 
9 
9 


9 
> 
9 


C 
NGO 
MIT 


C 
NGO 
MIT 


RAYN Ol 
RAYN 02 
RAYN 03 
RAYN 04 
RAYN O05 
RAYN O06 
RAYN O7 
RAYN 08 
RAYN O09 
RAYN 10 
RAYN 11 
RAYN 12 
RAYN 13 
RAYN 14 
RAYN 15 
RAYN 16 
RAYN 17 
RARAYN 18 
RAYN 19 
RAYN 20 
RAYN 21 
RAYN 22 
RAYN 23 
RAYN 24 
RAYN 25 
RAYN 26 
RAYN 27 
SURFCEO1 
SURFCEO2 
SURF CE03 
SURFCE0O4 
SURFCEO5 
SURFCE06 
SURFCEO7 
SURF CE0O8 
SURFCE09 
SURFCE10 
SURFCE11 
SURFCE12 
SURFCE13 
SURFCE14 
SURFCE15 
SURFCE16 
SURFCE17 
SURFCE18 
SURFCE19 
SURFCE20 
SURFCE21 


318 


DO 318 II=193 

YVW(TT) = Oo 

DO 318 L=194 

YVWCIT) = YVWCTI)+C(L)*EM(LoTT) 
DO 319 If=1s3 

E(II) = Oc 

DO 319 JJ=193 

E(TI) = ECIII+S(TIsJI)*#YVW( JJ) 
Oar 2 esl} se ECLQIIe se (EUS) AME 
FK = (E(2)*SIN(A)-E(3)*®COS(A))/CXKY 
RETURN 

END 


S$IBFTC MOVE 


102 
104 
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an Ww 


38 


SUBROUTINE MOVE (XsYoAsFK) 

DIMENSION S(393) 9EM( 493) 9E(3) sYVW(3) sCMAT( 100082) 9C(4&) oXLABL(12) 
COMMON S » EM 90 f » YVW » CMAT mG 
COMMON XLABL » D 5 Wi 9 CXY > IT » NGO 
COMMON AMM » ANN » MAX 9 NOT oN » MIT 
IF (MAX = 2) 10291025104 

FKBAR=FK 

MIT = 1 

DO 20 IT=1920 

DELA=FKBAR*®D 

AA=A4+DELA 

ABAR=A+e5*DELA 

DELX = D*COS(ABAR) 

DELY = D*SIN(ABAR) 

XX=X+DELX 

YY=Y+DELY 

GO TO (10196) MIT 

CALL SURFCE (XX9YYsAAsFKK) 

IF (CXY) 38938910 

FKBAR = 0¢5 * (FK + FKK) 

IF (IT — 18) 5293799 

FKKPP = FKBAR 

IF (MAX = 2) 79799 

IF (IT = 1) 2092009 

TF (ABS (FKKP-FKBAR) -— (0000009/D)) 696920 

FKKP = FKBAR 

IF (ABS (FKKPP = FKBAR) -— (02¢00009/D)) 18518917 
MIT = 3 

GO TO 38 

FKBAR = O0e5 * (FKBAR + FKKP) 

MIT = 2 

GO TO 39 

NGO = 1 = 

IF ((XX—100) #( (AMM—100)—-XX) 129293 

IF ((YY—-100)*( (ANN=100)-YY))29298 

NGO = 2 

X = XX 

V6 BAA 

A= AA 

FK = FKK 

RETURN 

END 
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SURFCE22 


SURFCE23 
SURFCE24 
SURFCE25 
SURFCE26 
SURFCE27 
SURFCE28 
SURFCE29 
SURF CE30 
SURFCE31 
SURFCE32 
SURFCE33 
MOVE 01 
MOVE 02 
MOVE 03 
MOVE 04 
MOVE 05 
MOVE 06 
MOVE 07 
MOVE 08 
MOVE 09 
MOVE 10 
MOVE 11 
MOVE 12 
MOVE 13 
MOVE 14 
MOVE 15 
MOVE 16 
MOVE 17 
MOVE 18 
MOVE 19 
MOVE 20 
MOVE 21 
MOVE 22 
MOVE 23 
MOVE 24 
MOVE 25 
MOVE 26 
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Derivations Relating QC with 2d, and 3C with 9 


Given: C= tanh | eTld 
” E ] 


Ox ayo MC 
Then: tanh} 2d} _ enc 
Cay - Ee 


efd _ tanh “1 fone 
Ge = 2 


d= ORs =o [iia iil Ave) 6 dia jl > Aue 
TI gt gr 


Y 
K 


Let: Teo cy TOG, eimgl KY Ss Piuyeae 
Then: @ = Gk? fin (@ & RG) sc tm (Go x"c)] 
Given: Cla) 6 (a) andardil— 21x59) 


However,it may also be condidered that: 
= iG) gael © S GOK, 99) 


Therefore (Kaplan, 1952, p.86): 


2d = d(d) . dC anddd _ d(a) . ac 
Tx 6a ox Ou) dot, os 
Sel — tlle re, - In(1-k"C)\9c +fin(1+k"C) - In(1-k"c)\oc 
Ox de Ox gs 
eG) a ie! ae y - (-k")\aC + fan(1+k"C) - In(1-k"c)\ac 
rp KC 1-k"CJQx ox 
2 -9d 1 
OX. “Soe ok! 
Similary: 
ac -gd. 1 iL 
QY oY &k' a 4 In(1+k"C) - In(1-k"C 
eC 


59 


APPENDIX G 
Method Used in Obtaining Most-Frequent 
Combinations of Deep-Water Height, Period, 

And Direction of Waves Capable of Striking Virginia Beach 

A rough estimate of the 15 most-frequent combinations of wave height, 
period, and direction of deep-water waves capable of striking Virginia Beach 
is gained by an analysis of five representative years of wave observations 
at the Chesapeake Lightship (fig. 1). Results of the analysis appear in 
table Gl and indicate that only six combinations of period and direction are 
involved in the 15 combinations. These six combinations are rough approxima- 
tions for input data, not only because of the crude methods involved in wave 
observation and recording, but also because it is assumed that one can pre- 
serve wave-direction angles for the waves observed at Chesapeake Lightship 
when searching for deep-water origin points that will yield rays capable of 
striking the target area. This is in error for some of the 4-second and 
all of the 6-second waves because they have already undergone refraction by 
the time they have reached Chesapeake Lightship; it is, therefore, invalid 
to select deep-water starting points for most waves by assuming that deep- 
water wave directions in water depths greater than 64 feet (Lightship water 
depth) are parallel to those observed at the Lightship. 

The 4 and 6-second waves that dominate the Lightship observations 
of table Gl are a product of the convention used here in interpreting the 
Ship's International Code. The convention is based upon an analysis of 
wave gage records that indicate a dominant 4-second period (table Ge) at 
a wave recording point off Cape Henry (fig. G1). (The Cape Henry wave 


gage was operated by the U. S. Navy between 1951 and 1956. It consisted 
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of a stepped-resistance wave staff; strip-chart recording were analyzed 
the significant wave method at 4-hour intervals.) Thus, code figure 2 
Code Table 17 of the Ship's International Code, which stands for waves of 
5-second period or less, was taken arbitrarily to represent waves of 


h4egecond period, while code figure 3 (5 to 7 seconds) was taken as 6 seconds. 
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Table Gl.-The fifteen most-frequently observed combinations of wave period, 
height, and direction (from 30°-150°T) as observed at Chesapeake 
Lightship during five representative years 


Combina- Period* Height Direction from No. of 
tion which wave front observa- 
(sec) (£t) travels (°T) tions** 
a 4 1 150 392 
2 4 3 60 286 
3 m 1 60 279 
4 4 1 90 265 
5 4 3 150 208 
6 4 3 90 166 
T 4 5 60 103 
8 6 3 60 100 
9 4 3 30 Qh 
10 6 3 90 93 
I. 4 1 30 We 
12 6 5 60 12 
13 6 al 90 50 
14 iF 6 60 ha 
i) in 5 30 39 


*Waves coded 2 {table 17, Ship's International Code) were considered as 
4-sec waves; those coded 3 were considered 6-sec waves. See explanation 
in text and table G2. 


**TOtal observations | of waves capable of striking Virginia Beach 
(traveling from 30° through 150° True) = 3016 
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FIGURE GI FORMER LOCATION OF U.S.NAVY WAVE GAGE OFF 
CAPE HENRY VIRGINIA, IN 20 FEET OF WATER 
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